#include <iostream>
#include <sstream>
#include <cstring>
#include <vector>
#include <queue>
#include <array>
#include <algorithm>
#include <tuple>
#include <unordered_map>
#include <fstream>
#include <cmath>
#include "getopt.h"
#include <sys/stat.h>
#include "config-fnt.h"
#include "graph.priv.h"
#include "stream.h"
#include "public/utils.h"

/* Options */
static const char* opts=":o:hv";
static const struct option opts_long[]=
{
	{"output", 1, nullptr, 'o'},
	{"help", 0, nullptr, 'h'},
	{"version", 0, nullptr, 'v'},
	{nullptr, 0, nullptr, 0}
};

/* Print version information and exit */
void version(const char* name) {
	std::cout<<name<<" (" PACKAGE_NAME_LONG ") " PACKAGE_VERSION "\n";
	std::cout<<"Copyright (C) 2015-2016 " PACKAGE_ORG "\n";
	std::cout<<"License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>\n\n";
	std::cout<<"Written by: " PACKAGE_AUTHORS "\n";
	exit(0);
}

/* Print usage information and exit */
void usage(const char* name, int ecode) {
	if(ecode) {
		std::cerr<<"Try '"<<name<<" --help' for more information.\n";
		exit(ecode);
	}
	std::cout<<"Usage: "<<name<<" [OPTION]... FNT1 FNT2\n";
	std::cout<<
		"Join FNT files.\n\n"
		"Options:\n"
		"  -o, --output=FILE [OPTIONAL] Save calculated distances between"
		"                       neurons to FILE. By default, save to standard "
		"                       output.\n\n"
		"  -h, --help        [OPTIONAL] Print this help, and do nothing else.\n"
		"  -v, --version     [OPTIONAL] Print version information and exit.\n\n"
		"Report bugs to <" PACKAGE_BUGREPORT ">\n"
		PACKAGE_NAME_LONG " home page: <" PACKAGE_URL ">\n";
	exit(ecode);
}

struct PointFloat {
	double x;
	double y;
	double z;
};
double dist(const PointFloat& p1, const PointFloat& p2) {
	double dx=p1.x-p2.x;
	double dy=p1.y-p2.y;
	double dz=p1.z-p2.z;
	return sqrt(dx*dx+dy*dy+dz*dz);
}
struct SerializedNeuron {
	struct Node {
		PointFloat pos;
		double distToParent;
		double lengthSum;
		size_t par;
		std::vector<size_t> branches;
		bool isbp;
		int16_t opt;
		double radius;
		Node() { }
		Node(const PointFloat& _pos, double d, double l, size_t p, const std::vector<size_t> br, bool b, int16_t o, double r=0): pos{_pos}, distToParent{d}, lengthSum{l},
			par{ p }, branches(br), isbp{ b }, opt{ o }, radius{ r } { }
	};
	std::vector<Node> nodes;
	SerializedNeuron(): nodes{} {
	}
	SerializedNeuron(size_t n): nodes{n} {
	}
	Node& node(size_t i) {
		return nodes[i];
	}
	const Node& node(size_t i) const {
		return nodes[i];
	}
};

std::vector<SerializedNeuron> serialize(Graph gr, double max_segm) {
	auto N=gr.trees().size();
	std::vector<SerializedNeuron> nrns(N);

	std::vector<bool> edgeVisited(gr.edges().size(), false);
	std::vector<size_t> vertNrn(gr.vertices().size(), SIZE_MAX);
	std::vector<size_t> vertIdx(gr.vertices().size(), SIZE_MAX);
	std::queue<Vertex> que{};

	for(size_t i=0; i<N; i++) {
		auto nrn=gr.trees()[i];
		auto soma=nrn.root();
		auto& nodes=nrns[i].nodes;
		if(soma.neighbors().size()==1) {
			vertNrn[soma.index()]=i;
			vertIdx[soma.index()]=SIZE_MAX;
			que.push(soma);
		} else {
			nodes.push_back({{soma.point().x(), soma.point().y(), soma.point().z()}, 0.0, 0.0, SIZE_MAX, {}, true, 1,soma.point().r()});
			vertNrn[soma.index()]=i;
			vertIdx[soma.index()]=0;
			que.push(soma);
		}
	}

	while(!que.empty()) {
		auto vert=que.front();
		que.pop();
		auto nrn=vertNrn[vert.index()];
		auto idx=vertIdx[vert.index()];
		auto& nodes=nrns[nrn].nodes;
		for(auto& p: vert.neighbors()) {
			auto edge=p.first;
			if(!edgeVisited[edge.index()]) {
				Vertex v1;
				auto parIdx=idx;
				if(p.second) {
					for(size_t j=edge.points().size(); j-->0;) {
						bool isbp=false;
						if(j==edge.points().size()-1) {
							isbp=false;
							if(parIdx!=SIZE_MAX)
								nodes[parIdx].branches.push_back(nodes.size());
						} else if(j==0) {
							isbp=true;
						}
						PointFloat pt{edge.points()[j].x(), edge.points()[j].y(), edge.points()[j].z()};
						if(j!=edge.points().size()-1) {
							auto pt0i=parIdx;
							auto l=dist(pt, nodes[pt0i].pos);
							if(l>max_segm) {
								int n=lrint(l/max_segm+.5);
								for(int i=1; i<n; i++) {
									PointFloat p{(nodes[pt0i].pos.x*(n-i)+pt.x*i)/n,(nodes[pt0i].pos.y*(n-i)+pt.y*i)/n,(nodes[pt0i].pos.z*(n-i)+pt.z*i)/n};
									nodes.push_back({p, 0.0, 0.0, parIdx, {}, false, edge.points()[j].m,edge.points()[j].r() });
									parIdx=nodes.size()-1;
								}
							}
						}
						nodes.push_back({pt, 0.0, 0.0, parIdx, {}, isbp, edge.points()[j].m,edge.points()[j].r() });
						parIdx=nodes.size()-1;
					}
					v1=edge.leftVertex();
				} else {
					for(size_t j=0; j<edge.points().size(); j++) {
						bool isbp=false;
						if(j==edge.points().size()-1) {
							isbp=true;
						} else if(j==0) {
							isbp=false;
							if(parIdx!=SIZE_MAX)
								nodes[parIdx].branches.push_back(nodes.size());
						}
						PointFloat pt{edge.points()[j].x(), edge.points()[j].y(), edge.points()[j].z()};
						if(j!=0) {
							auto pt0i=parIdx;
							auto l=dist(pt, nodes[pt0i].pos);
							if(l>max_segm) {
								int n=lrint(l/max_segm+.5);
								for(int k=1; k<n; k++) {
									PointFloat p{(nodes[pt0i].pos.x*(n-k)+pt.x*k)/n,(nodes[pt0i].pos.y*(n-k)+pt.y*k)/n,(nodes[pt0i].pos.z*(n-k)+pt.z*k)/n};
									nodes.push_back({p, 0.0, 0.0, parIdx, {}, false, edge.points()[j].m,edge.points()[j].r() });
									parIdx=nodes.size()-1;
								}
							}
						}
						nodes.push_back({{edge.points()[j].x(), edge.points()[j].y(), edge.points()[j].z()}, 0.0, 0.0, parIdx, {}, isbp, edge.points()[j].m,edge.points()[j].r() });
						parIdx=nodes.size()-1;
					}
					v1=edge.rightVertex();
				}
				edgeVisited[edge.index()]=true;
				if(vertNrn[v1.index()]==SIZE_MAX) {
					vertNrn[v1.index()]=nrn;
					vertIdx[v1.index()]=parIdx;
					que.push(v1);
				}
			}
		}
	}

	for(size_t i=0; i<N; i++) {
		auto& nodes=nrns[i].nodes;
		for(size_t j=nodes.size(); j-->0;) {
			auto par=nodes[j].par;
			if(par!=SIZE_MAX) {
				nodes[j].distToParent=dist(nodes[j].pos, nodes[par].pos);
			} else {
				nodes[j].distToParent=0;
			}

			double lengthSum=0;
			if(nodes[j].isbp) {
				for(auto b: nodes[j].branches) {
					lengthSum+=nodes[b].lengthSum+nodes[b].distToParent;
				}
			} else {
				lengthSum=nodes[j+1].lengthSum+nodes[j+1].distToParent;
			}
			nodes[j].lengthSum=lengthSum;
		}
	}
	return nrns;
}

struct Grid {
	size_t n0, n1;

	struct GridNode {
		double score;
		double sumMatch;
		double sumNomatch;
		size_t bti, btj;
		std::vector<std::pair<size_t, size_t>> bt;
		double dist;
		double scoreNoMatch0;
		double scoreNoMatch1;
		double mdist0;
		double mdist1;
	};
	std::vector<GridNode> nodes;

	Grid(size_t _n0, size_t _n1): n0{_n0}, n1{_n1}, nodes{n0*n1} {
	}
	GridNode& node(size_t i, size_t j) {
		return nodes[i*n1+j];
	}
};

void calcDist1(const PointFloat& a, const PointFloat& b, const PointFloat& c, const PointFloat& d, double l0, double l1, double* d0, double* d1) {
	auto abx=b.x-a.x;
	auto aby=b.y-a.y;
	auto abz=b.z-a.z;
	auto acx=c.x-a.x;
	auto acy=c.y-a.y;
	auto acz=c.z-a.z;
	auto ab2=abx*abx+aby*aby+abz*abz;
	auto abac=abx*acx+aby*acy+abz*acz;
	auto t0=abac/l0/l0;
	if(std::isnan(t0))
		t0=0;
	if(t0<0) t0=0;
	if(t0>1) t0=1;
	auto d02=t0*t0*l0*l0+ab2-2*t0*abac;
	if(d02<0) d02=0;
	*d0=sqrt(d02);
	auto bdx=d.x-b.x;
	auto bdy=d.y-b.y;
	auto bdz=d.z-b.z;
	auto abbd=abx*bdx+aby*bdy+abz*bdz;
	auto t1=-abbd/l1/l1;
	if(std::isnan(t1))
		t1=0;
	if(t1<0) t1=0;
	if(t1>1) t1=1;
	auto d12=t1*t1*l1*l1+ab2+2*t1*abbd;
	if(d12<0) d12=0;
	*d1=sqrt(d12);
}
struct NodePair {
	size_t i;
	size_t j;
};
struct NodePairScore {
	size_t i;
	size_t j;
	double s;
};
typedef std::tuple<int64_t, int16_t, double, double, double, double, int64_t> SwcNode;
typedef std::vector<SwcNode> SwcVector;
SwcVector calcDist(const SerializedNeuron& nrn0, const SerializedNeuron& nrn1, double* score, double* match, double* nomatch) {
	auto n0=nrn0.nodes.size();
	auto n1=nrn1.nodes.size();
	Grid grid{n0, n1};

	std::vector<NodePairScore> nodePairScores;
	std::vector<NodePair> nodePairs;
	std::vector<bool> nodeUsed0;
	std::vector<bool> nodeUsed1;

	for(size_t i=n0; i-->0;) {
		for(size_t j=n1; j-->0;) {
			auto& node=grid.node(i, j);
			node.dist=dist(nrn0.node(i).pos, nrn1.node(j).pos);
			auto paridx0=nrn0.node(i).par;
			if(paridx0==SIZE_MAX)
				paridx0=i;
			auto paridx1=nrn1.node(j).par;
			if(paridx1==SIZE_MAX)
				paridx1=j;
			calcDist1(nrn0.node(i).pos, nrn1.node(j).pos, nrn0.node(paridx0).pos, nrn1.node(paridx1).pos, nrn0.node(i).distToParent, nrn1.node(j).distToParent, &node.mdist0, &node.mdist1);
			double scoreNoMatch0=0;
			if(nrn0.node(i).isbp) {
				for(auto b: nrn0.node(i).branches)
					//scoreNoMatch0+=nrn0.node(b).distToParent*(node.dist+grid.node(b, j).dist)/2+grid.node(b, j).scoreNoMatch0;
					scoreNoMatch0+=/*nrn0.node(b).distToParent*grid.node(b, j).mdist0+*/grid.node(b, j).scoreNoMatch0;
			} else {
				//scoreNoMatch0+=nrn0.node(i+1).distToParent*(node.dist+grid.node(i+1, j).dist)/2+grid.node(i+1, j).scoreNoMatch0;
				scoreNoMatch0+=nrn0.node(i+1).distToParent*grid.node(i+1, j).mdist0+grid.node(i+1, j).scoreNoMatch0;
			}
			node.scoreNoMatch0=scoreNoMatch0;
			double scoreNoMatch1=0;
			if(nrn1.node(j).isbp) {
				for(auto b: nrn1.node(j).branches)
					//scoreNoMatch1+=nrn1.node(b).distToParent*(node.dist+grid.node(i, b).dist)/2+grid.node(i, b).scoreNoMatch1;
					scoreNoMatch1+=/*nrn1.node(b).distToParent*grid.node(i, b).mdist1+*/grid.node(i, b).scoreNoMatch1;
			} else {
				//scoreNoMatch1+=nrn1.node(j+1).distToParent*(node.dist+grid.node(i, j+1).dist)/2+grid.node(i, j+1).scoreNoMatch1;
				scoreNoMatch1+=nrn1.node(j+1).distToParent*grid.node(i, j+1).mdist1+grid.node(i, j+1).scoreNoMatch1;
			}
			node.scoreNoMatch1=scoreNoMatch1;

			node.score=INFINITY;

			if(!nrn0.node(i).isbp) {
				auto len=nrn0.node(i+1).distToParent;
				//auto mval=grid.node(i+1, j).score+len*(node.dist+grid.node(i+1,j).dist)/2;
				auto mval=grid.node(i+1, j).score+len*grid.node(i+1,j).mdist0;
				if(mval<node.score) {
					node.score=mval;
					//node.bt.resize(1);
					//node.bt[0].first=i+1;
					//node.bt[0].second=j;
					node.bti=i+1;
					node.btj=j;
					node.sumMatch=grid.node(i+1, j).sumMatch+len;
					node.sumNomatch=grid.node(i+1, j).sumNomatch;
				}
			}
			if(!nrn1.node(j).isbp) {
				auto len=nrn1.node(j+1).distToParent;
				//auto mval=grid.node(i, j+1).score+len*(node.dist+grid.node(i,j+1).dist)/2;
				auto mval=grid.node(i, j+1).score+len*grid.node(i,j+1).mdist1;
				if(mval<node.score) {
					node.score=mval;
					//node.bt.resize(1);
					//node.bt[0].first=i;
					//node.bt[0].second=j+1;
					node.bti=i;
					node.btj=j+1;
					node.sumMatch=grid.node(i, j+1).sumMatch+len;
					node.sumNomatch=grid.node(i, j+1).sumNomatch;
				}
				if(nrn0.node(i).isbp) {
					double dscore=INFINITY;
					size_t i1=SIZE_MAX;
					for(auto c: nrn0.node(i).branches) {
						double scored=grid.node(c, j).score-grid.node(c, j).scoreNoMatch0-grid.node(i, j).scoreNoMatch1;
						if(scored<dscore) {
							dscore=scored;
							i1=c;
						}
					}
					auto score=node.scoreNoMatch0+node.scoreNoMatch1+dscore;
					if(score<node.score) {
						node.score=score;
						//node.bt.resize(1);
						//node.bt[0].first=i1;
						//node.bt[0].second=j;
						node.bti=i1;
						node.btj=j;
						node.sumMatch=grid.node(i1, j).sumMatch;
						node.sumNomatch=grid.node(i1, j).sumNomatch+nrn0.node(i).lengthSum-nrn0.node(i1).lengthSum;
					}
				}
			} else {
				if(nrn0.node(i).isbp) {
					nodePairScores.clear();
					for(size_t i1p=0; i1p<nrn0.node(i).branches.size(); i1p++) {
						for(size_t j1p=0; j1p<nrn1.node(j).branches.size(); j1p++) {
							auto i1=nrn0.node(i).branches[i1p];
							auto j1=nrn1.node(j).branches[j1p];
							auto scored=grid.node(i1, j1).score-grid.node(i1, j).scoreNoMatch0-grid.node(i, j1).scoreNoMatch1;
							nodePairScores.push_back({i1p, j1p, scored});
						}
					}
					std::sort(nodePairScores.begin(), nodePairScores.end(), [](const NodePairScore& a, const NodePairScore& b) { return a.s<b.s; });

					double dscore=0;
					nodePairs.clear();
					nodeUsed0.assign(nrn0.node(i).branches.size(), false);
					nodeUsed1.assign(nrn1.node(j).branches.size(), false);
					for(auto& nps: nodePairScores) {
						if(nps.s>0)
							break;
						if(nodeUsed0[nps.i])
							continue;
						if(nodeUsed1[nps.j])
							continue;
						dscore+=nps.s;
						nodePairs.push_back({nps.i, nps.j});
						nodeUsed0[nps.i]=true;
						nodeUsed1[nps.j]=true;
					}

					if(dscore<=0) {
						node.score=node.scoreNoMatch0+node.scoreNoMatch1+dscore;
						node.bti=node.btj=SIZE_MAX;
						node.bt.resize(nodePairs.size());
						double sumMatch=0;
						double sumNomatch=nrn0.node(i).lengthSum+nrn1.node(j).lengthSum;
						for(size_t k=0; k<nodePairs.size(); k++) {
							auto i1=nrn0.node(i).branches[nodePairs[k].i];
							auto j1=nrn1.node(j).branches[nodePairs[k].j];
							node.bt[k].first=i1;
							node.bt[k].second=j1;
							sumMatch+=grid.node(i1, j1).sumMatch;
							sumNomatch+=grid.node(i1, j1).sumNomatch-nrn0.node(i1).lengthSum-nrn1.node(j1).lengthSum;
						}
						node.sumMatch=sumMatch;
						node.sumNomatch=sumNomatch;
					}
				} else {
					double dscore=INFINITY;
					size_t j1=SIZE_MAX;
					for(auto c: nrn1.node(j).branches) {
						auto scored=grid.node(i, c).score-grid.node(i, j).scoreNoMatch0-grid.node(i, c).scoreNoMatch1;
						if(scored<dscore) {
							dscore=scored;
							j1=c;
						}
					}
					auto score=node.scoreNoMatch0+node.scoreNoMatch1+dscore;
					if(dscore<0) {
						node.score=score;
						//node.bt.resize(1);
						//node.bt[0].first=i;
						//node.bt[0].second=j1;
						node.bti=i;
						node.btj=j1;
						node.sumMatch=grid.node(i, j1).sumMatch;
						node.sumNomatch=grid.node(i, j1).sumNomatch+nrn1.node(j).lengthSum-nrn1.node(j1).lengthSum;
					}
				}
			}
			if(node.scoreNoMatch0+node.scoreNoMatch1<node.score) {
				node.score=node.scoreNoMatch0+node.scoreNoMatch1;
				node.sumMatch=0;
				node.sumNomatch=nrn0.node(i).lengthSum+nrn1.node(j).lengthSum;
				node.bt.clear();
				node.bti=node.btj=SIZE_MAX;
			}
		}
	}

	struct TmpNode {
		int64_t id;
		int64_t par;
		size_t i;
		size_t j;
		bool samei;
		bool samej;
		int64_t xxx;
	};

	int64_t id=0;
	std::queue<TmpNode> que{};
	que.push(TmpNode{id++, -1, 0, 0, true, true, -1});
	std::vector<TmpNode> tmpnodes{};
	std::unordered_map<size_t, int64_t> i2node;
	std::unordered_map<size_t, int64_t> j2node;
	while(!que.empty()) {
		auto cur=que.front();
		que.pop();
		tmpnodes.push_back(cur);
		{
			auto it=i2node.find(cur.i);
			if(it==i2node.end())
				i2node[cur.i]=cur.id;
			it=j2node.find(cur.j);
			if(it==j2node.end())
				j2node[cur.j]=cur.id;
		}
		auto& node=grid.node(cur.i, cur.j);
		if(node.bti!=SIZE_MAX && node.btj!=SIZE_MAX) {
			que.push(TmpNode{id++, cur.id, node.bti, node.btj, true, true, -1});
		} else {
			for(auto& v: node.bt) {
				que.push(TmpNode{id++, cur.id, v.first, v.second, true, true, -1});
			}
		}
	}
	for(size_t i=0; i<nrn0.nodes.size(); i++) {
		auto it=i2node.find(i);
		if(it==i2node.end()) {
			i2node[i]=id;
			auto ip=i2node[nrn0.node(i).par];
			tmpnodes.push_back(TmpNode{id++, ip, i, tmpnodes[ip].j, true, true, -1});
		}
	}
	for(size_t j=0; j<nrn1.nodes.size(); j++) {
		auto it=j2node.find(j);
		if(it==j2node.end()) {
			j2node[j]=id;
			auto jp=j2node[nrn1.node(j).par];
			tmpnodes.push_back(TmpNode{id++, jp, tmpnodes[jp].i, j, true, true, -1});
		}
	}
	////////////

	for(size_t k=tmpnodes.size(); k-->0;) {
		auto kk=tmpnodes[k].par;
		if(kk>=0) {
			if(!tmpnodes[k].samei || tmpnodes[k].i!=tmpnodes[kk].i)
				tmpnodes[kk].samei=false;
			if(!tmpnodes[k].samej || tmpnodes[k].j!=tmpnodes[kk].j)
				tmpnodes[kk].samej=false;
		}
	}
	for(size_t k=0; k<tmpnodes.size(); k++) {
		if(tmpnodes[k].samei)
			tmpnodes[k].xxx=k;
		if(tmpnodes[k].samej)
			tmpnodes[k].xxx=k;
		auto kk=tmpnodes[k].par;
		if(kk>=0) {
			if(tmpnodes[kk].samei)
				tmpnodes[k].xxx=tmpnodes[kk].xxx;
			if(tmpnodes[kk].samej)
				tmpnodes[k].xxx=tmpnodes[kk].xxx;
		}
	}

	SwcVector swc(tmpnodes.size());
	for(size_t k=0; k<tmpnodes.size(); k++) {
		auto& n=tmpnodes[k];
		double x, y, z;
		int16_t opt=0;
		bool left=false, right=false;
		if(n.samei && n.samej) {
			auto kk=tmpnodes[k].par;
			auto& nn=tmpnodes[kk];
			if(nn.i==n.i)
				left=true;
			else if(nn.j==n.j)
				right=true;
		} else if(n.samei && !n.samej) {
			left=true;
		} else if(!n.samei && n.samej) {
			right=true;
		}
		float r= 0;
		if(left) {
			x=nrn1.node(n.j).pos.x;
			y=nrn1.node(n.j).pos.y;
			z=nrn1.node(n.j).pos.z;
			auto i=tmpnodes[n.xxx].i;
			auto j=tmpnodes[n.xxx].j;
			r = (nrn1.node(n.j).radius);
			x+=(nrn0.node(i).pos.x-nrn1.node(j).pos.x)/2;
			y+=(nrn0.node(i).pos.y-nrn1.node(j).pos.y)/2;
			z+=(nrn0.node(i).pos.z-nrn1.node(j).pos.z)/2;
			opt=7;
		} else if(right) {
			x=nrn0.node(n.i).pos.x;
			y=nrn0.node(n.i).pos.y;
			z=nrn0.node(n.i).pos.z;
			auto i=tmpnodes[n.xxx].i;
			auto j=tmpnodes[n.xxx].j;
			r = (nrn0.node(n.i).radius );
			x+=(nrn1.node(j).pos.x-nrn0.node(i).pos.x)/2;
			y+=(nrn1.node(j).pos.y-nrn0.node(i).pos.y)/2;
			z+=(nrn1.node(j).pos.z-nrn0.node(i).pos.z)/2;
			opt=7;
		} else {
			x=(nrn1.node(n.j).pos.x+nrn0.node(n.i).pos.x)/2;
			y=(nrn1.node(n.j).pos.y+nrn0.node(n.i).pos.y)/2;
			z=(nrn1.node(n.j).pos.z+nrn0.node(n.i).pos.z)/2;
			r = (nrn0.node(n.i).radius + nrn1.node(n.j).radius) / 2.0 +1;

			if(nrn1.node(n.j).opt==7)
				opt=6;
			if(nrn0.node(n.i).opt==7)
				opt=6;
		}
		swc[k]=SwcNode{n.id, opt, x, y, z, r, n.par};
	}
	std::get<1>(swc[0])=1;
	*score=grid.node(0, 0).score;
	*match=grid.node(0, 0).sumMatch;
	*nomatch=grid.node(0, 0).sumNomatch;
	/*
		for(auto& t: swc) {
			printMessage("asdf: ", std::get<0>(t));
		}
		*/
	return swc;
}

struct SimpleNode {
	double x, y;
	size_t par;
};

int main(int argc, char* argv[]) {
	const char* output=nullptr;
	double max_segm=INFINITY;

	// Parse options
	int opt;
	int prev_optind=optind;
	while((opt=getopt_long(argc, argv, opts, opts_long, nullptr))!=-1) {
		switch(opt) {
			case '?':
				std::cerr<<argv[0]<<": unrecognized option '";
				std::cerr<<argv[prev_optind]<<"'\n";
				usage(argv[0], -1);
			case ':':
				std::cerr<<argv[0]<<": option '"<<argv[prev_optind];
				std::cerr<<"' requires an argument\n";
				usage(argv[0], -1);
			case 'h':
				usage(argv[0], 0);
			case 'v':
				version(argv[0]);
			case 'o':
				output=optarg;
				break;
			default:
				std::cerr<<"Unexpected error.\n";
				usage(argv[0], -1);
		}
		prev_optind=optind;
	}


	if(optind+2!=argc) {
		std::cerr<<"No input FNT files specified.\n";
		usage(argv[0], -1);
	}

	if(!output) {
		std::cerr<<"No output file given\n";
	}

	Graph gr1=GraphPriv::alloc();
	auto gp1=GraphPriv::get(gr1);
	auto headers=gp1->loadFnt(argv[optind]);

	Graph gr2=GraphPriv::alloc();
	auto gp2=GraphPriv::get(gr2);
	headers=gp2->loadFnt(argv[optind+1]);

	auto neurons1=serialize(gr1, max_segm);
	auto neurons2=serialize(gr2, max_segm);
	if(neurons1.size()!=1 || neurons2.size()!=1) {
		std::cerr<<"please give fnt files, each containing one neuron.\n";
		return -1;
	}



	double dist, match, nomatch;
	auto swc=calcDist(neurons1[0], neurons2[0], &dist, &match, &nomatch);
	GraphPriv::free(gr1);
	GraphPriv::free(gr2);

	try {
		Graph gr=GraphPriv::alloc();
		auto gp=GraphPriv::get(gr);
		gp->fromSwc(swc);
		gp->saveFnt(output, headers);
		GraphPriv::free(gr);
	} catch(std::exception& e) {
		std::cerr<<"Error while saving FNT file: "<<e.what();
		return -1;
	}

	return 0;
}

